Problem 1 : Tomato

Remember the tomato data set generated by Pepe; we first looked at this when we were working on Chapter 9. 35 accessions for seven species were grown in sun and shade.

Assess whether there is evidence for total length (“totleng”) response to shade and whether this response varies by species. Consider whether including accession (“acs”), and/or shelf (“shelf”) using adaptive priors improves the model fit.

Bonus: would it be better to consider shade by accession interactions instead of shade x species?

setwd("~/R /R_club/R_Club_2016/Stacey.Harmer_RethinkingHomework/Assigment_2016-12-12")
tom <- read.csv("TomatoR2CSHL.csv")
head(tom)
##   shelf flat col row    acs trt days   date   hyp int1 int2 int3 int4
## 1     Z    1   B   1 LA2580   H   28 5/5/08 19.46 2.37 1.59 1.87 0.51
## 2     Z    1   C   1 LA1305   H   28 5/5/08 31.28 3.34 0.01 9.19 1.62
## 3     Z    1   D   1 LA1973   H   28 5/5/08 56.65 8.43 2.39 6.70 3.69
## 4     Z    1   E   1 LA2748   H   28 5/5/08 35.18 0.56 0.00 1.60 0.61
## 5     Z    1   F   1 LA2931   H   28 5/5/08 35.32 0.82 0.02 1.49 0.46
## 6     Z    1   G   1 LA1317   H   28 5/5/08 28.74 1.07 6.69 5.72 4.76
##   intleng totleng petleng leafleng leafwid leafnum ndvi      lat      lon
## 1    6.34   25.80   15.78    30.53   34.44       5  111  -9.5167 -78.0083
## 2   14.16   45.44   12.36    22.93   13.99       4  120 -13.3833 -75.3583
## 3   21.21   77.86   13.05    46.71   43.78       5  110 -16.2333 -71.7000
## 4    2.77   37.95    8.08    26.82   33.28       5  105 -20.4833 -69.9833
## 5    2.79   38.11    7.68    22.40   23.61       5  106 -20.9167 -69.0667
## 6   18.24   46.98   23.66    42.35   42.35       5  132 -13.4167 -73.8417
##    alt         species who
## 1  740    S. pennellii Dan
## 2 3360   S. peruvianum Dan
## 3 2585   S. peruvianum Dan
## 4 1020     S. chilense Dan
## 5 2460     S. chilense Dan
## 6 2000 S. chmielewskii Dan
str(tom)
## 'data.frame':    1008 obs. of  25 variables:
##  $ shelf   : Factor w/ 6 levels "U","V","W","X",..: 6 6 6 6 6 6 6 6 6 6 ...
##  $ flat    : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ col     : Factor w/ 9 levels "A","B","C","D",..: 2 3 4 5 6 7 9 2 3 4 ...
##  $ row     : int  1 1 1 1 1 1 1 2 2 2 ...
##  $ acs     : Factor w/ 36 levels "LA1028","LA1305",..: 23 2 19 26 30 5 14 32 7 34 ...
##  $ trt     : Factor w/ 2 levels "H","L": 1 1 1 1 1 1 1 1 1 1 ...
##  $ days    : int  28 28 28 28 28 28 28 28 28 28 ...
##  $ date    : Factor w/ 2 levels "5/5/08","5/6/08": 1 1 1 1 1 1 1 1 1 1 ...
##  $ hyp     : num  19.5 31.3 56.6 35.2 35.3 ...
##  $ int1    : num  2.37 3.34 8.43 0.56 0.82 1.07 2.85 2.08 5.43 4.08 ...
##  $ int2    : num  1.59 0.01 2.39 0 0.02 6.69 0.41 0.53 0.81 3.26 ...
##  $ int3    : num  1.87 9.19 6.7 1.6 1.49 5.72 3.79 1.9 3.63 3.49 ...
##  $ int4    : num  0.51 1.62 3.69 0.61 0.46 4.76 3.25 NA 2.66 3.02 ...
##  $ intleng : num  6.34 14.16 21.21 2.77 2.79 ...
##  $ totleng : num  25.8 45.4 77.9 38 38.1 ...
##  $ petleng : num  15.78 12.36 13.05 8.08 7.68 ...
##  $ leafleng: num  30.5 22.9 46.7 26.8 22.4 ...
##  $ leafwid : num  34.4 14 43.8 33.3 23.6 ...
##  $ leafnum : int  5 4 5 5 5 5 5 4 5 5 ...
##  $ ndvi    : int  111 120 110 105 106 132 118 112 107 123 ...
##  $ lat     : num  -9.52 -13.38 -16.23 -20.48 -20.92 ...
##  $ lon     : num  -78 -75.4 -71.7 -70 -69.1 ...
##  $ alt     : int  740 3360 2585 1020 2460 2000 2920 480 75 3540 ...
##  $ species : Factor w/ 5 levels "S. chilense",..: 4 5 5 1 1 2 3 4 5 5 ...
##  $ who     : Factor w/ 2 levels "Dan","Pepe": 1 1 1 1 1 1 1 1 1 1 ...
summary(tom)
##  shelf        flat            col           row            acs     
##  U:161   Min.   : 1.00   G      :133   Min.   :1.00   LA1954 : 40  
##  V:174   1st Qu.: 9.00   H      :127   1st Qu.:2.00   LA2695 : 39  
##  W:178   Median :17.00   F      :125   Median :3.00   LA1361 : 37  
##  X:174   Mean   :17.89   C      :117   Mean   :2.56   LA2167 : 37  
##  Y:125   3rd Qu.:28.00   D      :117   3rd Qu.:4.00   LA2773 : 37  
##  Z:196   Max.   :36.00   E      :107   Max.   :4.00   LA1474 : 36  
##                          (Other):282                  (Other):782  
##  trt          days           date          hyp             int1      
##  H:495   Min.   :28.00   5/5/08:716   Min.   : 6.17   Min.   : 0.00  
##  L:513   1st Qu.:28.00   5/6/08:292   1st Qu.:26.81   1st Qu.: 1.74  
##          Median :28.00                Median :32.02   Median : 3.59  
##          Mean   :28.29                Mean   :33.36   Mean   : 4.71  
##          3rd Qu.:29.00                3rd Qu.:38.56   3rd Qu.: 6.46  
##          Max.   :29.00                Max.   :74.60   Max.   :39.01  
##                                                       NA's   :1      
##       int2             int3             int4           intleng      
##  Min.   : 0.000   Min.   : 0.010   Min.   : 0.030   Min.   : 0.000  
##  1st Qu.: 1.060   1st Qu.: 2.975   1st Qu.: 2.163   1st Qu.: 9.637  
##  Median : 3.120   Median : 5.625   Median : 3.995   Median :17.255  
##  Mean   : 4.287   Mean   : 6.794   Mean   : 5.102   Mean   :20.340  
##  3rd Qu.: 6.320   3rd Qu.: 9.367   3rd Qu.: 7.018   3rd Qu.:28.145  
##  Max.   :28.980   Max.   :27.760   Max.   :23.280   Max.   :92.420  
##  NA's   :1        NA's   :4        NA's   :102                      
##     totleng          petleng         leafleng        leafwid     
##  Min.   : 13.59   Min.   : 1.53   Min.   : 9.74   Min.   : 8.29  
##  1st Qu.: 39.25   1st Qu.:11.20   1st Qu.:27.43   1st Qu.:29.48  
##  Median : 50.98   Median :15.13   Median :34.59   Median :39.62  
##  Mean   : 53.70   Mean   :15.92   Mean   :35.54   Mean   :39.29  
##  3rd Qu.: 64.76   3rd Qu.:20.48   3rd Qu.:42.98   3rd Qu.:47.75  
##  Max.   :129.43   Max.   :44.44   Max.   :95.19   Max.   :90.27  
##                   NA's   :2       NA's   :1       NA's   :1      
##     leafnum           ndvi          lat               lon        
##  Min.   :3.000   Min.   :100   Min.   :-25.400   Min.   :-78.52  
##  1st Qu.:5.000   1st Qu.:108   1st Qu.:-16.607   1st Qu.:-75.92  
##  Median :5.000   Median :115   Median :-14.152   Median :-73.63  
##  Mean   :5.063   Mean   :118   Mean   :-14.490   Mean   :-73.71  
##  3rd Qu.:6.000   3rd Qu.:128   3rd Qu.:-12.450   3rd Qu.:-71.70  
##  Max.   :8.000   Max.   :137   Max.   : -5.767   Max.   :-68.07  
##  NA's   :1                                                       
##       alt                  species      who     
##  Min.   :   0   S. chilense    :207   Dan :402  
##  1st Qu.:1020   S. chmielewskii:226   Pepe:606  
##  Median :2240   S. habrochaites:226             
##  Mean   :2035   S. pennellii   :132             
##  3rd Qu.:3110   S. peruvianum  :217             
##  Max.   :3540                                   
## 
library(rethinking)
## Loading required package: rstan
## Loading required package: ggplot2
## Loading required package: StanHeaders
## rstan (Version 2.12.1, packaged: 2016-09-11 13:07:50 UTC, GitRev: 85f7a56811da)
## For execution on a local, multicore CPU with excess RAM we recommend calling
## rstan_options(auto_write = TRUE)
## options(mc.cores = parallel::detectCores())
## Loading required package: parallel
## rethinking (Version 1.59)
# get map2stan up and ready
library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

Make dummy variable, so that sun = 0, shade = 1

tom$shade <- ifelse(tom$trt == "L", 1, 0)
head(tom, n=6)
##   shelf flat col row    acs trt days   date   hyp int1 int2 int3 int4
## 1     Z    1   B   1 LA2580   H   28 5/5/08 19.46 2.37 1.59 1.87 0.51
## 2     Z    1   C   1 LA1305   H   28 5/5/08 31.28 3.34 0.01 9.19 1.62
## 3     Z    1   D   1 LA1973   H   28 5/5/08 56.65 8.43 2.39 6.70 3.69
## 4     Z    1   E   1 LA2748   H   28 5/5/08 35.18 0.56 0.00 1.60 0.61
## 5     Z    1   F   1 LA2931   H   28 5/5/08 35.32 0.82 0.02 1.49 0.46
## 6     Z    1   G   1 LA1317   H   28 5/5/08 28.74 1.07 6.69 5.72 4.76
##   intleng totleng petleng leafleng leafwid leafnum ndvi      lat      lon
## 1    6.34   25.80   15.78    30.53   34.44       5  111  -9.5167 -78.0083
## 2   14.16   45.44   12.36    22.93   13.99       4  120 -13.3833 -75.3583
## 3   21.21   77.86   13.05    46.71   43.78       5  110 -16.2333 -71.7000
## 4    2.77   37.95    8.08    26.82   33.28       5  105 -20.4833 -69.9833
## 5    2.79   38.11    7.68    22.40   23.61       5  106 -20.9167 -69.0667
## 6   18.24   46.98   23.66    42.35   42.35       5  132 -13.4167 -73.8417
##    alt         species who shade
## 1  740    S. pennellii Dan     0
## 2 3360   S. peruvianum Dan     0
## 3 2585   S. peruvianum Dan     0
## 4 1020     S. chilense Dan     0
## 5 2460     S. chilense Dan     0
## 6 2000 S. chmielewskii Dan     0
tail(tom, n=20)
##      shelf flat col row    acs trt days   date   hyp int1 int2 int3 int4
## 989      X   36   A   2 LA2167   H   28 5/5/08 30.00 1.71 3.00 4.16 4.46
## 990      X   36   B   2 LA3788   H   28 5/5/08 19.70 0.62 0.74 0.90 0.43
## 991      X   36   C   2 LA1336   H   28 5/5/08 35.19 3.43 2.64 5.93 5.27
## 992      X   36   D   2 LA3795   H   28 5/5/08 26.17 2.06 1.99 0.87 5.65
## 993      X   36   E   2 LA2773   H   28 5/5/08 13.02 1.50 1.13 2.79 2.78
## 994      X   36   F   2 LA3788   H   28 5/5/08 23.39 5.47 7.44 9.95 4.74
## 995      X   36   H   2 LA1361   H   28 5/5/08 26.34 1.75 2.46 5.42 3.88
## 996      X   36   A   3 LA1969   H   28 5/5/08 26.36 0.16 0.61 0.30 0.13
## 997      X   36   C   3 LA1474   H   28 5/5/08 34.95 2.87 2.58 5.27 4.28
## 998      X   36   E   3 LA2884   H   28 5/5/08 26.32 1.35 0.44 3.46 2.11
## 999      X   36   F   3 LA1306   H   28 5/5/08 27.06 1.86 5.96 4.97 3.00
## 1000     X   36   G   3 LA2663   H   28 5/5/08 26.03 1.42 4.27 5.35 3.97
## 1001     X   36   H   3 LA1363   H   28 5/5/08 30.61 3.87 4.44 4.49 2.90
## 1002     X   36   I   3 LA2167   H   28 5/5/08 25.04 4.30 2.16 2.18   NA
## 1003     X   36   A   4 LA1474   H   28 5/5/08 29.38 0.12 0.44 0.44 0.21
## 1004     X   36   D   4 LA1969   H   28 5/5/08 37.00 1.01 0.72 0.70 1.20
## 1005     X   36   F   4 LA1316   H   28 5/5/08 26.59 2.00 5.17 6.06 3.01
## 1006     X   36   G   4 LA2695   H   28 5/5/08 32.11 1.26 2.29 3.69 4.55
## 1007     X   36   H   4 LA2695   H   28 5/5/08 29.00 3.06 7.37 6.14 4.12
## 1008     X   36   I   4 LA2204   H   28 5/5/08 27.62 2.03 4.99 3.75   NA
##      intleng totleng petleng leafleng leafwid leafnum ndvi      lat
## 989    13.33   43.33   12.42    49.90   51.83       6  124  -7.1800
## 990     2.69   22.39    8.03    23.08   24.92       4  112 -16.1383
## 991    17.27   52.46   16.86    33.21   43.20       6  107 -16.2094
## 992    10.57   36.74      NA    21.63   16.77       4  123 -10.1500
## 993     8.20   21.22   10.34    26.83   25.09       5  115 -18.2667
## 994    27.60   50.99   21.52    45.31   50.90       4  112 -16.1383
## 995    13.51   39.85   16.08    22.49   26.75       5  111  -9.5500
## 996     1.20   27.56    3.25    20.38   19.19       5  124 -17.5333
## 997    15.00   49.95   16.87    40.00   55.34       6  113 -16.6072
## 998     7.36   33.68    9.06    37.37   39.86       5  108 -23.6833
## 999    15.79   42.85   13.33    31.02   36.17       4  130 -12.9500
## 1000   15.01   41.04   16.12    40.15   48.18       5  137 -13.7333
## 1001   15.70   46.31   15.32    31.76   40.30       4  121 -10.1417
## 1002    8.64   33.68   13.94    16.96   22.22       5  124  -7.1800
## 1003    1.21   30.59    4.35    14.56   14.19       4  113 -16.6072
## 1004    3.63   40.63    6.42    18.67   23.17       5  124 -17.5333
## 1005   16.24   42.83   16.48    30.80   30.52       4  135 -13.3925
## 1006   11.79   43.90   23.77    42.93   50.53       5  128 -14.1525
## 1007   20.69   49.69   15.08    38.75   47.90       5  128 -14.1525
## 1008   10.77   38.39   10.25    22.97   26.65       4  128  -5.7667
##           lon  alt         species  who shade
## 989  -78.5200 1340 S. habrochaites Pepe     0
## 990  -73.6428  480    S. pennellii Pepe     0
## 991  -73.6222   75   S. peruvianum Pepe     0
## 992  -77.3500 3540   S. peruvianum Pepe     0
## 993  -69.5833 3260     S. chilense Pepe     0
## 994  -73.6428  480    S. pennellii Pepe     0
## 995  -77.8917 1340 S. habrochaites Pepe     0
## 996  -70.0333 3015     S. chilense Pepe     0
## 997  -72.6964   27   S. peruvianum Pepe     0
## 998  -68.0667 2390     S. chilense Pepe     0
## 999  -74.0167 3340 S. chmielewskii Pepe     0
## 1000 -71.9833 3110 S. chmielewskii Pepe     0
## 1001 -77.3833 3140 S. habrochaites Pepe     0
## 1002 -78.5200 1340 S. habrochaites Pepe     0
## 1003 -72.6964   27   S. peruvianum Pepe     0
## 1004 -70.0333 3015     S. chilense Pepe     0
## 1005 -73.9150 3100 S. chmielewskii Pepe     0
## 1006 -71.3886 3480 S. chmielewskii Pepe     0
## 1007 -71.3886 3480 S. chmielewskii Pepe     0
## 1008 -77.8667 2160 S. habrochaites Pepe     0

I want to treat species as a cluster variable. Do I need to give each species a number? Try without doing that.

First write model only taking treatment and species into account After first attempt, learned that I need to trim this dataset. Keep accession (acs) and shelf and totleng and treatment and speices but ditch the rest

tom.trim <- tom[, c(1,5, 15, 24,26)]
   
head(tom.trim)
##   shelf    acs totleng         species shade
## 1     Z LA2580   25.80    S. pennellii     0
## 2     Z LA1305   45.44   S. peruvianum     0
## 3     Z LA1973   77.86   S. peruvianum     0
## 4     Z LA2748   37.95     S. chilense     0
## 5     Z LA2931   38.11     S. chilense     0
## 6     Z LA1317   46.98 S. chmielewskii     0
str(tom.trim)
## 'data.frame':    1008 obs. of  5 variables:
##  $ shelf  : Factor w/ 6 levels "U","V","W","X",..: 6 6 6 6 6 6 6 6 6 6 ...
##  $ acs    : Factor w/ 36 levels "LA1028","LA1305",..: 23 2 19 26 30 5 14 32 7 34 ...
##  $ totleng: num  25.8 45.4 77.9 38 38.1 ...
##  $ species: Factor w/ 5 levels "S. chilense",..: 4 5 5 1 1 2 3 4 5 5 ...
##  $ shade  : num  0 0 0 0 0 0 0 0 0 0 ...
tom.trim$sp_index <- coerce_index(tom.trim$species)
mTom.1 <- map2stan(
  alist(
    totleng ~ dnorm(mu, sigma),
    mu <- a_species[sp_index] + b_species[sp_index]*shade,
    a_species[sp_index] ~ dnorm(0, 10),
    b_species[sp_index] ~ dnorm(0,10),
    sigma ~ dnorm(0, 10)), data=tom.trim, iter=5000, warmup=2000, chains =2)
## Warning in FUN(X[[i]], ...): data with name shelf is not numeric and not
## used
## Warning in FUN(X[[i]], ...): data with name acs is not numeric and not used
## Warning in FUN(X[[i]], ...): data with name species is not numeric and not
## used
## 
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 1).
## 
## Chain 1, Iteration:    1 / 5000 [  0%]  (Warmup)
## Chain 1, Iteration:  500 / 5000 [ 10%]  (Warmup)
## Chain 1, Iteration: 1000 / 5000 [ 20%]  (Warmup)
## Chain 1, Iteration: 1500 / 5000 [ 30%]  (Warmup)
## Chain 1, Iteration: 2000 / 5000 [ 40%]  (Warmup)
## Chain 1, Iteration: 2001 / 5000 [ 40%]  (Sampling)
## Chain 1, Iteration: 2500 / 5000 [ 50%]  (Sampling)
## Chain 1, Iteration: 3000 / 5000 [ 60%]  (Sampling)
## Chain 1, Iteration: 3500 / 5000 [ 70%]  (Sampling)
## Chain 1, Iteration: 4000 / 5000 [ 80%]  (Sampling)
## Chain 1, Iteration: 4500 / 5000 [ 90%]  (Sampling)
## Chain 1, Iteration: 5000 / 5000 [100%]  (Sampling)
##  Elapsed Time: 1.75539 seconds (Warm-up)
##                1.36139 seconds (Sampling)
##                3.11678 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 2).
## 
## Chain 2, Iteration:    1 / 5000 [  0%]  (Warmup)
## Chain 2, Iteration:  500 / 5000 [ 10%]  (Warmup)
## Chain 2, Iteration: 1000 / 5000 [ 20%]  (Warmup)
## Chain 2, Iteration: 1500 / 5000 [ 30%]  (Warmup)
## Chain 2, Iteration: 2000 / 5000 [ 40%]  (Warmup)
## Chain 2, Iteration: 2001 / 5000 [ 40%]  (Sampling)
## Chain 2, Iteration: 2500 / 5000 [ 50%]  (Sampling)
## Chain 2, Iteration: 3000 / 5000 [ 60%]  (Sampling)
## Chain 2, Iteration: 3500 / 5000 [ 70%]  (Sampling)
## Chain 2, Iteration: 4000 / 5000 [ 80%]  (Sampling)
## Chain 2, Iteration: 4500 / 5000 [ 90%]  (Sampling)
## Chain 2, Iteration: 5000 / 5000 [100%]  (Sampling)
##  Elapsed Time: 1.55028 seconds (Warm-up)
##                1.34519 seconds (Sampling)
##                2.89547 seconds (Total)
## Warning in FUN(X[[i]], ...): data with name shelf is not numeric and not
## used
## Warning in FUN(X[[i]], ...): data with name acs is not numeric and not used
## Warning in FUN(X[[i]], ...): data with name species is not numeric and not
## used
## 
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 1).
## WARNING: No variance estimation is
##          performed for num_warmup < 20
## 
## 
## Chain 1, Iteration: 1 / 1 [100%]  (Sampling)
##  Elapsed Time: 4e-06 seconds (Warm-up)
##                0.000142 seconds (Sampling)
##                0.000146 seconds (Total)
## Computing WAIC
## Constructing posterior predictions
## [ 600 / 6000 ]
[ 1200 / 6000 ]
[ 1800 / 6000 ]
[ 2400 / 6000 ]
[ 3000 / 6000 ]
[ 3600 / 6000 ]
[ 4200 / 6000 ]
[ 4800 / 6000 ]
[ 5400 / 6000 ]
[ 6000 / 6000 ]
plot(mTom.1)
precis(mTom.1, depth=2)
##               Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a_species[1] 36.88   1.40      34.66      39.10  4317    1
## a_species[2] 43.83   1.33      41.77      45.98  4190    1
## a_species[3] 44.34   1.30      42.26      46.42  4272    1
## a_species[4] 29.07   1.82      26.18      32.04  4082    1
## a_species[5] 48.20   1.36      46.02      50.33  3977    1
## b_species[1] 28.03   1.91      25.07      31.15  4164    1
## b_species[2] 17.10   1.83      14.10      19.86  4203    1
## b_species[3] 20.40   1.88      17.29      23.32  4192    1
## b_species[4] 26.62   2.42      22.74      30.41  4197    1
## b_species[5] 24.66   1.89      21.79      27.77  4522    1
## sigma        15.03   0.31      14.53      15.52  6000    1
plot(precis(mTom.1, depth = 2))

Consider whether including accession (“acs”), and/or shelf (“shelf”) using adaptive priors improves the model fit.

# OK, include accession and then shelf; first make them an index

tom.trim$acs_index <- coerce_index(tom.trim$acs)
tom.trim$shelf_index <- coerce_index(tom.trim$shelf)
head(tom.trim)
##   shelf    acs totleng         species shade sp_index acs_index
## 1     Z LA2580   25.80    S. pennellii     0        4        23
## 2     Z LA1305   45.44   S. peruvianum     0        5         2
## 3     Z LA1973   77.86   S. peruvianum     0        5        19
## 4     Z LA2748   37.95     S. chilense     0        1        26
## 5     Z LA2931   38.11     S. chilense     0        1        30
## 6     Z LA1317   46.98 S. chmielewskii     0        2         5
##   shelf_index
## 1           6
## 2           6
## 3           6
## 4           6
## 5           6
## 6           6
# include accession, using adaptive prior
mTom.2 <- map2stan(
  alist(
    totleng ~ dnorm(mu, sigma),
    mu <- a_species[sp_index] + a_acs[acs_index] + b_species[sp_index]*shade,
    a_species[sp_index] ~ dnorm(0, 10),
    b_species[sp_index] ~ dnorm(0,10),
    a_acs[acs_index] ~ dnorm(0, sigma_acs),
    sigma_acs ~ dnorm(0, 10),
    sigma ~ dnorm(0, 10)), data=tom.trim, iter=5000, warmup=2000, chains =2)
## Warning in FUN(X[[i]], ...): data with name shelf is not numeric and not
## used
## Warning in FUN(X[[i]], ...): data with name acs is not numeric and not used
## Warning in FUN(X[[i]], ...): data with name species is not numeric and not
## used
## 
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 1).
## 
## Chain 1, Iteration:    1 / 5000 [  0%]  (Warmup)
## Chain 1, Iteration:  500 / 5000 [ 10%]  (Warmup)
## Chain 1, Iteration: 1000 / 5000 [ 20%]  (Warmup)
## Chain 1, Iteration: 1500 / 5000 [ 30%]  (Warmup)
## Chain 1, Iteration: 2000 / 5000 [ 40%]  (Warmup)
## Chain 1, Iteration: 2001 / 5000 [ 40%]  (Sampling)
## Chain 1, Iteration: 2500 / 5000 [ 50%]  (Sampling)
## Chain 1, Iteration: 3000 / 5000 [ 60%]  (Sampling)
## Chain 1, Iteration: 3500 / 5000 [ 70%]  (Sampling)
## Chain 1, Iteration: 4000 / 5000 [ 80%]  (Sampling)
## Chain 1, Iteration: 4500 / 5000 [ 90%]  (Sampling)
## Chain 1, Iteration: 5000 / 5000 [100%]  (Sampling)
##  Elapsed Time: 6.35734 seconds (Warm-up)
##                5.93577 seconds (Sampling)
##                12.2931 seconds (Total)
## The following numerical problems occured the indicated number of times after warmup on chain 1
##                                                                                 count
## Exception thrown at line 21: normal_log: Scale parameter is 0, but must be > 0!     2
## When a numerical problem occurs, the Hamiltonian proposal gets rejected.
## See http://mc-stan.org/misc/warnings.html#exception-hamiltonian-proposal-rejected
## If the number in the 'count' column is small, do not ask about this message on stan-users.
## 
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 2).
## 
## Chain 2, Iteration:    1 / 5000 [  0%]  (Warmup)
## Chain 2, Iteration:  500 / 5000 [ 10%]  (Warmup)
## Chain 2, Iteration: 1000 / 5000 [ 20%]  (Warmup)
## Chain 2, Iteration: 1500 / 5000 [ 30%]  (Warmup)
## Chain 2, Iteration: 2000 / 5000 [ 40%]  (Warmup)
## Chain 2, Iteration: 2001 / 5000 [ 40%]  (Sampling)
## Chain 2, Iteration: 2500 / 5000 [ 50%]  (Sampling)
## Chain 2, Iteration: 3000 / 5000 [ 60%]  (Sampling)
## Chain 2, Iteration: 3500 / 5000 [ 70%]  (Sampling)
## Chain 2, Iteration: 4000 / 5000 [ 80%]  (Sampling)
## Chain 2, Iteration: 4500 / 5000 [ 90%]  (Sampling)
## Chain 2, Iteration: 5000 / 5000 [100%]  (Sampling)
##  Elapsed Time: 5.7553 seconds (Warm-up)
##                4.98947 seconds (Sampling)
##                10.7448 seconds (Total)
## The following numerical problems occured the indicated number of times after warmup on chain 2
##                                                                                 count
## Exception thrown at line 21: normal_log: Scale parameter is 0, but must be > 0!     2
## When a numerical problem occurs, the Hamiltonian proposal gets rejected.
## See http://mc-stan.org/misc/warnings.html#exception-hamiltonian-proposal-rejected
## If the number in the 'count' column is small, do not ask about this message on stan-users.
## Warning in FUN(X[[i]], ...): data with name shelf is not numeric and not
## used
## Warning in FUN(X[[i]], ...): data with name acs is not numeric and not used
## Warning in FUN(X[[i]], ...): data with name species is not numeric and not
## used
## 
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 1).
## WARNING: No variance estimation is
##          performed for num_warmup < 20
## 
## 
## Chain 1, Iteration: 1 / 1 [100%]  (Sampling)
##  Elapsed Time: 3e-06 seconds (Warm-up)
##                0.000202 seconds (Sampling)
##                0.000205 seconds (Total)
## Computing WAIC
## Constructing posterior predictions
## [ 600 / 6000 ]
[ 1200 / 6000 ]
[ 1800 / 6000 ]
[ 2400 / 6000 ]
[ 3000 / 6000 ]
[ 3600 / 6000 ]
[ 4200 / 6000 ]
[ 4800 / 6000 ]
[ 5400 / 6000 ]
[ 6000 / 6000 ]
plot(mTom.2)
## Waiting to draw page 2 of 4

## Waiting to draw page 3 of 4

## Waiting to draw page 4 of 4

# neff not great for the a_acs values
precis(mTom.2, depth=2)
##               Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a_species[1] 32.80   3.80      27.24      39.24  1146    1
## a_species[2] 39.05   3.93      33.09      45.30  1015    1
## a_species[3] 39.53   3.66      33.79      45.32   879    1
## a_species[4] 26.24   3.72      20.41      32.04  1214    1
## a_species[5] 42.67   3.83      37.18      48.86  1090    1
## b_species[1] 27.67   1.82      24.96      30.72  4889    1
## b_species[2] 16.83   1.74      14.07      19.56  4961    1
## b_species[3] 19.98   1.73      17.23      22.75  4424    1
## b_species[4] 26.98   2.26      23.70      30.86  4112    1
## b_species[5] 24.35   1.76      21.56      27.17  5099    1
## a_acs[1]     -1.06   4.51      -7.83       6.32  1404    1
## a_acs[2]     -6.38   4.24     -13.13       0.13  1412    1
## a_acs[3]      1.10   4.25      -5.96       7.24  1179    1
## a_acs[4]      7.99   4.35       0.78      14.34  1187    1
## a_acs[5]     18.75   4.40      11.49      25.19  1134    1
## a_acs[6]      7.29   4.36       0.63      14.21  1204    1
## a_acs[7]     14.95   4.34       8.18      21.63  1300    1
## a_acs[8]     10.98   4.08       4.88      17.55  1173    1
## a_acs[9]      8.03   4.11       1.50      14.44   982    1
## a_acs[10]     8.24   4.52       1.28      15.59  1648    1
## a_acs[11]     7.17   4.22       0.51      13.60  1269    1
## a_acs[12]    -5.72   4.09     -12.36       0.49  1410    1
## a_acs[13]     0.17   4.18      -6.72       6.35  1508    1
## a_acs[14]     3.65   4.29      -2.98      10.65  1297    1
## a_acs[15]    -0.65   5.43      -9.12       8.15  3129    1
## a_acs[16]     6.23   4.19      -0.80      12.32  1247    1
## a_acs[17]     9.35   4.20       2.61      15.76  1380    1
## a_acs[18]    -8.07   4.08     -14.19      -1.46  1527    1
## a_acs[19]    20.27   4.46      13.20      27.16  1169    1
## a_acs[20]    -1.86   3.99      -7.98       4.61  1197    1
## a_acs[21]     5.47   4.02      -0.46      12.28  1250    1
## a_acs[22]    10.57   5.15       2.18      18.52  1880    1
## a_acs[23]     2.00   4.23      -4.99       8.25  1587    1
## a_acs[24]     0.36   4.25      -5.87       7.44  1189    1
## a_acs[25]    -0.46   4.21      -7.30       5.91  1174    1
## a_acs[26]     4.07   4.38      -2.33      11.49  1481    1
## a_acs[27]     6.38   4.15      -0.44      12.61  1334    1
## a_acs[28]    -0.33   4.16      -6.99       6.19  1468    1
## a_acs[29]     5.56   4.31      -1.43      12.14  1447    1
## a_acs[30]    11.45   4.28       4.86      18.31  1354    1
## a_acs[31]    15.15   4.27       8.57      22.11  1213    1
## a_acs[32]     4.85   4.18      -1.96      11.16  1407    1
## a_acs[33]     0.94   4.11      -5.09       7.89  1501    1
## a_acs[34]    -4.13   4.22     -10.99       2.14  1337    1
## a_acs[35]    -1.17   4.31      -8.31       5.15  1418    1
## a_acs[36]    -3.99   4.79     -11.94       3.36  2180    1
## sigma_acs     9.25   1.79       6.62      11.81   934    1
## sigma        13.34   0.30      12.86      13.81  5213    1
par(mfrow = c(1,1))

plot(precis(mTom.2, depth = 2))

# looks a bit worse than the simpler model

# include shelf, using adaptive prior
mTom.3 <- map2stan(
  alist(
    totleng ~ dnorm(mu, sigma),
    mu <- a_species[sp_index] + a_shelf[shelf_index] + b_species[sp_index]*shade,
    a_species[sp_index] ~ dnorm(0, 10),
    b_species[sp_index] ~ dnorm(0,10),
    a_shelf[shelf_index] ~ dnorm(0, sigma_acs),
    sigma_acs ~ dnorm(0, 10),
    sigma ~ dnorm(0, 10)), data=tom.trim, iter=5000, warmup=2000, chains =2)
## Warning in FUN(X[[i]], ...): data with name shelf is not numeric and not
## used
## Warning in FUN(X[[i]], ...): data with name acs is not numeric and not used
## Warning in FUN(X[[i]], ...): data with name species is not numeric and not
## used
## 
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 1).
## 
## Chain 1, Iteration:    1 / 5000 [  0%]  (Warmup)
## Chain 1, Iteration:  500 / 5000 [ 10%]  (Warmup)
## Chain 1, Iteration: 1000 / 5000 [ 20%]  (Warmup)
## Chain 1, Iteration: 1500 / 5000 [ 30%]  (Warmup)
## Chain 1, Iteration: 2000 / 5000 [ 40%]  (Warmup)
## Chain 1, Iteration: 2001 / 5000 [ 40%]  (Sampling)
## Chain 1, Iteration: 2500 / 5000 [ 50%]  (Sampling)
## Chain 1, Iteration: 3000 / 5000 [ 60%]  (Sampling)
## Chain 1, Iteration: 3500 / 5000 [ 70%]  (Sampling)
## Chain 1, Iteration: 4000 / 5000 [ 80%]  (Sampling)
## Chain 1, Iteration: 4500 / 5000 [ 90%]  (Sampling)
## Chain 1, Iteration: 5000 / 5000 [100%]  (Sampling)
##  Elapsed Time: 7.10618 seconds (Warm-up)
##                5.47194 seconds (Sampling)
##                12.5781 seconds (Total)
## The following numerical problems occured the indicated number of times after warmup on chain 1
##                                                                                 count
## Exception thrown at line 27: normal_log: Scale parameter is 0, but must be > 0!     4
## When a numerical problem occurs, the Hamiltonian proposal gets rejected.
## See http://mc-stan.org/misc/warnings.html#exception-hamiltonian-proposal-rejected
## If the number in the 'count' column is small, do not ask about this message on stan-users.
## 
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 2).
## 
## Chain 2, Iteration:    1 / 5000 [  0%]  (Warmup)
## Chain 2, Iteration:  500 / 5000 [ 10%]  (Warmup)
## Chain 2, Iteration: 1000 / 5000 [ 20%]  (Warmup)
## Chain 2, Iteration: 1500 / 5000 [ 30%]  (Warmup)
## Chain 2, Iteration: 2000 / 5000 [ 40%]  (Warmup)
## Chain 2, Iteration: 2001 / 5000 [ 40%]  (Sampling)
## Chain 2, Iteration: 2500 / 5000 [ 50%]  (Sampling)
## Chain 2, Iteration: 3000 / 5000 [ 60%]  (Sampling)
## Chain 2, Iteration: 3500 / 5000 [ 70%]  (Sampling)
## Chain 2, Iteration: 4000 / 5000 [ 80%]  (Sampling)
## Chain 2, Iteration: 4500 / 5000 [ 90%]  (Sampling)
## Chain 2, Iteration: 5000 / 5000 [100%]  (Sampling)
##  Elapsed Time: 6.25795 seconds (Warm-up)
##                7.16284 seconds (Sampling)
##                13.4208 seconds (Total)
## Warning in FUN(X[[i]], ...): data with name shelf is not numeric and not
## used
## Warning in FUN(X[[i]], ...): data with name acs is not numeric and not used
## Warning in FUN(X[[i]], ...): data with name species is not numeric and not
## used
## 
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 1).
## WARNING: No variance estimation is
##          performed for num_warmup < 20
## 
## 
## Chain 1, Iteration: 1 / 1 [100%]  (Sampling)
##  Elapsed Time: 4e-06 seconds (Warm-up)
##                0.000207 seconds (Sampling)
##                0.000211 seconds (Total)
## Computing WAIC
## Constructing posterior predictions
## [ 600 / 6000 ]
[ 1200 / 6000 ]
[ 1800 / 6000 ]
[ 2400 / 6000 ]
[ 3000 / 6000 ]
[ 3600 / 6000 ]
[ 4200 / 6000 ]
[ 4800 / 6000 ]
[ 5400 / 6000 ]
[ 6000 / 6000 ]
plot(mTom.3)
## Waiting to draw page 2 of 2

precis(mTom.3, depth=2) # again, n-eff not great
##               Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a_species[1]  2.93   4.48      -4.15      10.19  1043    1
## a_species[2]  9.76   4.47       2.74      17.16  1062    1
## a_species[3] 10.36   4.48       3.38      17.58  1053    1
## a_species[4] -5.56   4.53     -12.73       1.69  1124    1
## a_species[5] 14.16   4.49       6.91      21.15  1077    1
## b_species[1]  8.21   4.64       0.94      15.45  1501    1
## b_species[2] -2.67   4.65     -10.12       4.68  1478    1
## b_species[3]  1.25   4.65      -6.36       8.28  1562    1
## b_species[4]  7.06   4.76      -0.50      14.77  1592    1
## b_species[5]  4.75   4.65      -2.89      11.81  1526    1
## a_shelf[1]   50.76   6.20      40.92      60.61   877    1
## a_shelf[2]   57.94   6.18      48.20      67.82   881    1
## a_shelf[3]   53.60   6.19      43.52      63.28   890    1
## a_shelf[4]   32.21   4.43      25.10      39.28  1027    1
## a_shelf[5]   32.95   4.47      25.71      40.05  1045    1
## a_shelf[6]   37.36   4.41      30.31      44.40  1037    1
## sigma_acs    30.15   4.70      22.75      37.39  2170    1
## sigma        14.82   0.33      14.29      15.33  3532    1
par(mfrow = c(1,1))

plot(precis(mTom.3, depth = 2))

# looks like shlef has a bigger effect than species

now include both acs and shelf in the model

mTom.4 <- map2stan(
  alist(
    totleng ~ dnorm(mu, sigma),
    mu <- a_species[sp_index] + a_shelf[shelf_index] + a_acs[acs_index] + b_species[sp_index]*shade,
    a_species[sp_index] ~ dnorm(0, 10),
    b_species[sp_index] ~ dnorm(0,10),
    a_acs[acs_index] ~ dnorm(0, sigma_acs),
    a_shelf[shelf_index] ~ dnorm(0, sigma_acs),
    sigma_acs ~ dnorm(0, 10),
    sigma_acs ~ dnorm(0, 10),
    sigma ~ dnorm(0, 10)), data=tom.trim, iter=5000, warmup=2000, chains =2)
## Warning in FUN(X[[i]], ...): data with name shelf is not numeric and not
## used
## Warning in FUN(X[[i]], ...): data with name acs is not numeric and not used
## Warning in FUN(X[[i]], ...): data with name species is not numeric and not
## used
## 
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 1).
## 
## Chain 1, Iteration:    1 / 5000 [  0%]  (Warmup)
## Chain 1, Iteration:  500 / 5000 [ 10%]  (Warmup)
## Chain 1, Iteration: 1000 / 5000 [ 20%]  (Warmup)
## Chain 1, Iteration: 1500 / 5000 [ 30%]  (Warmup)
## Chain 1, Iteration: 2000 / 5000 [ 40%]  (Warmup)
## Chain 1, Iteration: 2001 / 5000 [ 40%]  (Sampling)
## Chain 1, Iteration: 2500 / 5000 [ 50%]  (Sampling)
## Chain 1, Iteration: 3000 / 5000 [ 60%]  (Sampling)
## Chain 1, Iteration: 3500 / 5000 [ 70%]  (Sampling)
## Chain 1, Iteration: 4000 / 5000 [ 80%]  (Sampling)
## Chain 1, Iteration: 4500 / 5000 [ 90%]  (Sampling)
## Chain 1, Iteration: 5000 / 5000 [100%]  (Sampling)
##  Elapsed Time: 11.0538 seconds (Warm-up)
##                15.1305 seconds (Sampling)
##                26.1844 seconds (Total)
## The following numerical problems occured the indicated number of times after warmup on chain 1
##                                                                                 count
## Exception thrown at line 25: normal_log: Scale parameter is 0, but must be > 0!     2
## When a numerical problem occurs, the Hamiltonian proposal gets rejected.
## See http://mc-stan.org/misc/warnings.html#exception-hamiltonian-proposal-rejected
## If the number in the 'count' column is small, do not ask about this message on stan-users.
## 
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 2).
## 
## Chain 2, Iteration:    1 / 5000 [  0%]  (Warmup)
## Chain 2, Iteration:  500 / 5000 [ 10%]  (Warmup)
## Chain 2, Iteration: 1000 / 5000 [ 20%]  (Warmup)
## Chain 2, Iteration: 1500 / 5000 [ 30%]  (Warmup)
## Chain 2, Iteration: 2000 / 5000 [ 40%]  (Warmup)
## Chain 2, Iteration: 2001 / 5000 [ 40%]  (Sampling)
## Chain 2, Iteration: 2500 / 5000 [ 50%]  (Sampling)
## Chain 2, Iteration: 3000 / 5000 [ 60%]  (Sampling)
## Chain 2, Iteration: 3500 / 5000 [ 70%]  (Sampling)
## Chain 2, Iteration: 4000 / 5000 [ 80%]  (Sampling)
## Chain 2, Iteration: 4500 / 5000 [ 90%]  (Sampling)
## Chain 2, Iteration: 5000 / 5000 [100%]  (Sampling)
##  Elapsed Time: 11.2639 seconds (Warm-up)
##                17.1731 seconds (Sampling)
##                28.4371 seconds (Total)
## The following numerical problems occured the indicated number of times after warmup on chain 2
##                                                                                 count
## Exception thrown at line 25: normal_log: Scale parameter is 0, but must be > 0!     2
## When a numerical problem occurs, the Hamiltonian proposal gets rejected.
## See http://mc-stan.org/misc/warnings.html#exception-hamiltonian-proposal-rejected
## If the number in the 'count' column is small, do not ask about this message on stan-users.
## Warning in FUN(X[[i]], ...): data with name shelf is not numeric and not
## used
## Warning in FUN(X[[i]], ...): data with name acs is not numeric and not used
## Warning in FUN(X[[i]], ...): data with name species is not numeric and not
## used
## 
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 1).
## WARNING: No variance estimation is
##          performed for num_warmup < 20
## 
## 
## Chain 1, Iteration: 1 / 1 [100%]  (Sampling)
##  Elapsed Time: 3e-06 seconds (Warm-up)
##                0.000225 seconds (Sampling)
##                0.000228 seconds (Total)
## Computing WAIC
## Constructing posterior predictions
## [ 600 / 6000 ]
[ 1200 / 6000 ]
[ 1800 / 6000 ]
[ 2400 / 6000 ]
[ 3000 / 6000 ]
[ 3600 / 6000 ]
[ 4200 / 6000 ]
[ 4800 / 6000 ]
[ 5400 / 6000 ]
[ 6000 / 6000 ]
plot(mTom.4) 
## Waiting to draw page 2 of 4

## Waiting to draw page 3 of 4

## Waiting to draw page 4 of 4

precis(mTom.4, depth=2) # n-eff better now
##               Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a_species[1] 15.01   5.92       5.41      24.13  2072    1
## a_species[2] 20.45   6.20      10.69      30.66  2132    1
## a_species[3] 20.98   6.14      11.09      30.52  1935    1
## a_species[4]  8.17   5.83      -1.08      17.53  2032    1
## a_species[5] 23.84   6.26      13.71      33.39  1913    1
## b_species[1] 15.38   4.55       7.76      22.17  2472    1
## b_species[2]  4.68   4.59      -2.52      12.14  2302    1
## b_species[3]  8.31   4.54       0.79      15.22  2404    1
## b_species[4] 15.21   4.73       8.00      23.01  2574    1
## b_species[5] 12.03   4.57       4.92      19.37  2483    1
## a_acs[1]     -2.61   5.33     -10.42       6.60  2990    1
## a_acs[2]     -7.15   5.05     -15.07       0.95  2617    1
## a_acs[3]      1.29   5.10      -6.54       9.48  2889    1
## a_acs[4]      8.39   5.06       0.13      16.26  2726    1
## a_acs[5]     19.46   5.12      11.08      27.45  2695    1
## a_acs[6]      7.64   5.11      -0.49      15.51  2762    1
## a_acs[7]     15.79   5.07       7.46      23.59  2369    1
## a_acs[8]     11.77   4.84       3.97      19.14  2743    1
## a_acs[9]      8.30   4.82       0.83      15.79  2631    1
## a_acs[10]     6.10   5.19      -2.48      14.06  3349    1
## a_acs[11]     8.13   4.98       0.38      16.31  2378    1
## a_acs[12]    -6.64   4.95     -14.66       1.06  2824    1
## a_acs[13]    -1.53   4.98      -9.88       6.03  3194    1
## a_acs[14]     3.41   5.01      -4.76      11.06  2738    1
## a_acs[15]    -3.95   6.46     -14.26       6.09  3994    1
## a_acs[16]     7.30   4.97      -0.39      15.51  2491    1
## a_acs[17]     9.30   4.79       1.80      16.90  2924    1
## a_acs[18]    -9.38   4.88     -16.56      -1.09  2993    1
## a_acs[19]    21.28   5.07      12.88      29.02  2373    1
## a_acs[20]    -0.87   4.80      -8.29       6.75  2628    1
## a_acs[21]     6.19   4.83      -1.37      13.88  2671    1
## a_acs[22]     7.83   5.91      -1.84      16.77  3677    1
## a_acs[23]     0.52   5.07      -7.81       7.99  3286    1
## a_acs[24]     0.56   5.12      -7.41       8.75  2885    1
## a_acs[25]     0.05   5.02      -7.60       8.18  2799    1
## a_acs[26]     3.64   5.08      -4.39      11.71  2732    1
## a_acs[27]     6.24   4.80      -1.26      13.99  2868    1
## a_acs[28]    -1.18   4.83      -8.68       6.58  2848    1
## a_acs[29]     5.27   4.91      -2.86      12.62  2776    1
## a_acs[30]    11.05   4.84       3.81      19.31  2697    1
## a_acs[31]    15.93   4.97       8.01      23.62  2774    1
## a_acs[32]     4.49   4.97      -3.43      12.44  3175    1
## a_acs[33]     0.17   4.94      -7.89       7.76  3098    1
## a_acs[34]    -4.03   4.97     -11.87       3.84  2565    1
## a_acs[35]    -1.33   5.13      -8.88       7.27  2284    1
## a_acs[36]    -4.60   5.58     -13.60       4.18  3471    1
## a_shelf[1]   26.82   6.35      16.60      36.78  1620    1
## a_shelf[2]   34.70   6.36      25.08      45.28  1613    1
## a_shelf[3]   30.37   6.33      20.27      40.36  1593    1
## a_shelf[4]   15.80   4.50       8.67      22.97  1891    1
## a_shelf[5]   17.34   4.55       9.99      24.54  1913    1
## a_shelf[6]   21.75   4.48      14.45      28.82  1867    1
## sigma_acs    13.06   2.09       9.87      16.38  1918    1
## sigma        13.04   0.29      12.60      13.50  5486    1
par(mfrow = c(1,1))

plot(precis(mTom.4, depth = 2))

# 

compare these 4 models

compare(mTom.1, mTom.2, mTom.3, mTom.4)
##          WAIC pWAIC dWAIC weight    SE   dSE
## mTom.4 8083.7  44.3   0.0      1 56.24    NA
## mTom.2 8124.1  39.2  40.4      0 55.24 14.20
## mTom.3 8310.2  14.5 226.5      0 54.42 31.29
## mTom.1 8334.7  10.2 251.0      0 53.84 32.95

As I suspected based on the n_eff values, model 4 is much preferred.

Problem 2: Smoking deaths among doctors

In 1961 Doll and Hill sent out a questionnaire to all men on the British Medical Register inquiring about their smoking habits. Almost 70% of such men replied. Death certificates were obtained for medical practitioners and causes of death were assigned on the basis of these certificates. The breslow data set contains the person-years of observations and deaths from coronary artery disease accumulated during the first ten years of the study.

Analyse this data set to determine the posterior probability that smoking increases death by coronary artery disease, that age increases death by coronary artery disease, and that there is an interaction between age and smoking.

You can load the data set and learn about the columns using the commands below

data(“breslow”,package = “boot”) help(“breslow”,package =“boot”) breslow You can think of “person-years” as the number of observations

Note: You almost certainly have the boot package on your computer, but if you do not have the boot package on your computer then you will need to install.packages(“boot”)

Note: do NOT do library(boot). This will make the logit function from boot over-ride the one from rethinking Note: You will probably need to do:

breslow\(n <- as.integer(breslow\)n)

before you can analyze the problem 2 data

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.

Problem 3: Cane Sugar

This data comes from an experiment to measure disease resistance in different varieties of sugar cane.

Is there evidence of differences in disease resistance in the different varieties? Does including an adaptive prior for Block improve the model fit?

You can get the data and learn about it with:

data(“cane”,package=“boot”) help(“cane”,package=“boot”) head(cane) summary(cane)